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Abstract 

The concept of isogeometric analysis, whereby the parametric ftmc- 
tions that are used to describe CAD geometry are also used to approx- 
imate the unknown fields in a numerical discretisation, has progressed 
rapidly in recent years. This paper advances the field further by outhn- 
ing an isogeometric Boundary Element Method (IGABEM) that only re- 
quires a representation of the geometry of the domain for analysis, fitting 
neatly with the boundary representation provided completely by CAD. 
The method circumvents the requirement to generate a boundary mesh 
representing a significant step in reducing the gap between engineering 
design and analysis. The current paper focuses on implementation details 
of 2D IGABEM for elastostatic analysis with particular attention paid 
towards the differences over conventional boundary element implementa- 
tions. Examples of Matlab® code are given whenever possible to aid 
understanding of the techniques used. 

1. Introduction 

Isogeometric analysis (IGA) is a subject which is receiving a great deal of 
attention amongst the computational mechanics community since it has the 
potential to have a profound effect on the current engineering design and anal- 
ysis process. The concept has the capability of leading to large steps forward 
in efficiency since effectively, the process of "meshing" is either eliminated or 
greatly suppressed. To understand why this is such a revolutionary develop- 
ment, it is necessary to take an abstract view of the current engineering design 
and analysis process and evaluate the critical points where bottlenecks occur 
which subsequently lead to delays in engineering projects. 
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Models are created in Computer Aided Design (CAD) software allowing 
designers to materialise ideas into computational objects that can range from 
simple geometries to highly detailed and complex engineering prototypes. The 
wide array of CAD packages available, along with the ever-increasing geometri- 
cal modelling capabilities, offers designers the ability to create realistic models 
of complex components. Once a model is complete in CAD, it must be trans- 
formed into a form suitable for analysis, with the most time-consuming step 
taken up in creating a suitable "analysis-ready" model which forms a discreti- 
sation of the domain (or boundary). This step in many cases requires human 
intervention by a specialist to ensure that the discretisation is of sufficient qual- 
ity to give accurate results in future simulations. But what is most important 
to note, is that the relative portion of time taken to create an analysis-ready 
design is approximately 80% of the total design and analysis process, thereby 
dominating the entire design process. 

Once a suitable discretisation has been made, analysis can be carried out 
using a suitable numerical method such as the Finite Element Method (FEM), 
Finite Difference Method (EDM) or the Boundary Element Method (BEM). 
Analysis itself represents a relatively small portion of the design process but 
importantly, we find that often the original design is affected by the results ob- 
tained from analysis. In this manner, design and analysis are tightly connected 
through an iterative procedure - a concept which is mirrored throughout the 
engineering community. 

Recently, an answer to the problems created by the mismatch between design 
and analysis was proposed through the concept of Isogeometric Analysis. The 
concept was initially proposed by Hughes et al. [l9| and since this seminal work, 
a book has been published entirely on the subject [l^l ■ Rather than using con- 
ventional piecewise polynomial shape functions to discretise both the geometry 
and unknown fields, IGA proposes to use the parametric functions used by CAD 
as an approximation for both fields, most commonly taking the form of Non- 
Uniform Rational B-Splines (NURBS). In this way, the isoparametric concept is 
maintained but more significantly, the geometry of the problem is preserved ex- 
actly. In addition, since many of the algorithms implemented in CAD packages 
can also be used for numerical analysis, redundant computations are eliminated 
allowing analysis to be carried out with greatly reduced pre-processing. 

A great deal of research has been focussed on IGA in recent years with 
implementations in areas such as patient-specific modelling 0, XFEM 



shells and many others e.g.[l|,[3|,[23|i[l3|i[23|'|2li- In the majority of these 
methods NURBS are used for discretisation, but the inability of the functions to 
produce "watertight" geometries and allow local refinement have shown major 
shortcomings. Perhaps one of the most significant developments to overcome the 
deficiencies of NURBS is the introduction of T-splines Q and later PHT-splines 
[2^ that produce watertight geometries and can be locally refined, benefiting 
both the design and numerical analysis communities. In particular, from an 
analysis standpoint, the use of such functions is essential for efficient algorithms 
to exploit adaptivity. More recently IGA has been applied within a BEM context 
for elastostatic analysis [2^ where particular benefits are realised due to the 
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requirement for only a surface representation of the geometry. The present 
paper builds on this work, where emphasis is given towards the implementation 
details of a 2D isogeometric BEM. The organisation of the paper is as follows: 
first, an outline of B-splines and NURBS which form the underlying technology 
of isogeometric analysis is given, a review of the conventional BEM is given 
and finally, the implementation details of isogeometric BEM are illustrated by 
building-up an example problem from an inital CAD model to the final BEM 
system of equations. 



2. Geometrical modelling 

The key concept of isogeometric analysis is bringing the fields of design 
and analysis together into a unified framework through the use of parametric 
functions that are predominant in CAD. Therefore we concentrate in the current 
section on describing such functions which most commonly take the form of 
NURBS. In much of the recent literature on isogeometric analysis 13| lojjsl 



(and indeed, much literature in the past ^245,^261), extensive details are given 
on the construction of B-splines and NURBS and therefore we only give a brief 
description of these functions ensuring that relevant notation is defined to aid 
the reader in later sections. 

Since the present paper is focussed on implementation details, heavy use is 
made of the algorithms stated in [24] which cover details ranging from evaluating 
NURBS basis functions to refinement algorithms such as knot insertion and 
order elevation. The reader is advised to consult this reference since the basics 
of parametric functions for geometric modelling are outlined clearly. 



2.1. B-splines and NURBS 

The first concept which must be grasped when using B-splines or NURBS 
is that both functions are parametric in that the equations which describe the 
curves (or surfaces) are completely defined by a number of independent param- 
eters. In the current context we will use ^ as the independent variable which 
describes a B-spline or NURBS and denote it as a coordinate in the parameter 
space. To understand the meaning of the variable ^, an example B-spline is 
shown in Fig[T] which illustrates the basic but important concepts of the func- 
tions that are used to decribe CAD geometry. We note the following: 

• The curve requires a set of control points to be defined which may or 
may not lie on the curve 

• The curve requires the definition of a knot vector, defined as a non- 
decreasing sequence of coordinates in the parameter space, which in this 
case is given by S = {0, 0, 0, 1, 2, 3, 4, 4, 4} 

• In this instance the curve is constructed from an open knot vector which 
results in a spline that is interpolatory at the beginning and end of the 
curve. 
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Figure 1: Example B-spline with associated control points. 



• If knot values are repeated, then it is found that the order of continuity 
decreases at that point in the spline. 

The knot vector is a concept which is often unfamilar to those working in the 
field of numerical analysis but should be considered carefully, since it has a large 
influence on the resulting spline. The most important aspect of the knot vector 
is the relative difference between the components and for this reason the values 
can be scaled if required. That is, the knot vector S = {0, 0, 0, 1, 2, 3, 4, 4, 4} with 
^ e [0, 4) is equivalent to S = {0, 0, 0, 0.1, 0.2, 0.3, 0.4, 0.4, 0.4} with ^ e [0, 0.4). 
In later sections we will see that the process of applying refinement in IGABEM 
has a considerable effect on the knot vector. Knot vectors are ubiqutous in 
the fields of geometrical modelling and isogeometric analysis and therefore the 
reader should become accustomed to their use. 

We now introduce some more formal definitions which are required for suc- 
cint notation in later sections. Denoting the dimensionality of the problem as 
R'' {d = 2, 3), the following three items are required to define fully a B-spline: 

• The curve degree, p, e.g. linear {p = 1), quadratic {p ~ 2)... 

• A set of of n control points Pa E M.'^ , I < a < n 

• A knot vector S = {^i, ^2, 6n+p-hi} 

Figures [5] to |3] illustrate linear, quadratic and cubic B-splines with their 
associated control points. The knot vectors associated with each curve are 



{0, 0, 1, 2, 3, 3}, {0, 0, 0, 1, 2, 3, 3, 3} and {0, 0, 0, 0, 1, 2, 3, 3, 3, 3} respectively. Each 
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Figure 4: Cubic B-spline, p = 3 
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of these are found to be open knot vectors which denotes that they contain p + 1 
repeated components at the beginning and end of the knot vector. The fact that 
the curve is interpolatory at the beginning and end is a direct consequence of 
this. In ah future examples used in this paper, it can be assumed that open knot 
vectors are used. Repeated knot components may also occur at points which 
are not located at the extremes with a consequence that the continuity of the 
curve is reduced at that point. 

2.1.1. B-spline and NURBS basis functions 

The previous section served as an overview of B-splines and some of the ter- 
minlogy associated with their construction. But for a B-spline to be completely 
defined, some attention must be paid towards their associated basis functions. 
The idea of interpolating a discrete number of points mirrors the technology seen 
in conventional FEMs and BEMs but with a distinct difference - the curve is 
not required to exhibit the Kronecker delta property at the interpolated points. 
The consequences of this are that when interpolating fields such as displacement, 
the value obtained at nodal points does not represent any real displacement but 
rather a coefficient used for interpolation. This is similar to meshless meth- 
ods where techniques such as Lagrange multipliers [8] and the penalty method 
approach Q are employed to impose essential boundary conditions. 

Turning our attention towards B-spline basis functions, we introduce the fol- 
lowing expression which fully describes a B-spline in terms of its basis functions 
and control points. This is written as 

n 
a=l 

where C(^) is a vector denoting the Cartesian coordinates of the location de- 
scribed by the parametric coordinate ^, and Na.p{£,) denotes the set of B-spline 
basis functions of degree p at ^. The basis functions are defined as 

if U<£.< L+i 
otherwise 




and for p = 1,2, 3. 



sa+p t,a 
, Ca+p+1 ~ S, 



-Na+l,p-liO- (3) 
?a+p+l — t;a+l 



Fig. [S] illustrates the basis functions corresponding to the B-splinc shown in 
Fig. [1] here the interpolatory nature of the first and last basis functions is evi- 
dent. What should be noted is that Eqns. ([2]) and (O are recursive in nature and, 
in their current form, considerably more expensive than conventional polynomial 
basis function expressions. However, there exist several efficient computational 
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Figure 5: B-spline basis functions for curve shown in Fig. [4] and knot vector H = 
{0,0,0,0,1,2,3,3,3,3} 



algorithms for their evaluation such as the Cox-de-Boor algorithm j24j . |26| and, 
more recently, the extraction operator [27]. B-spline derivatives are also re- 
quired for numerical analysis, but the algorithms for determining their values 
are standard in CAD literature, with details given in [Appendix A.l[ 

In CAD surface modelling packages, NURBS represent the dominant tool 
used to describe curves and surfaces where, in fact, they are found to be a 
superset of B-splines and only differ from their counterparts by the use of an 
additional coordinate often referred to as a 'weighting'. In some interpretations 
NURBS are seen as a 'projection' of B-splines from a higher dimensional space 
(see [ij]), and it can be shown that some attractive properties emerge. In par- 
ticular, NURBS are able to reproduce circular arcs, spheres and conic sections 
exactly (cf. B-splines which only approximate such shapes) and this is achieved 
through the appropriate choice of weightings. Each control point is associ- 
ated with a weighting Wa leading to a set of NURBS basis functions denoted by 
Ra,p{£,)- The curve is then interpolated as 

n 

C(0-E^-.p(OPa (4) 

a=l 

with the NURBS basis functions given by 
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Figure 6: Definition of problem domain witii source and field points. 



where Na^p can be found from Eqns ^ and In the case that the weights are 
all set to unity (i.e. = 1 Va), the basis functions given by ([5]) reduce to the 
B-spline basis functions of ([3]) . Expressions for the derivaties of NURBS basis 
functions can also be found, and arc detailed in [Appendix A.2| Since B-spUnes 
are a subset of NURBS, for the sake of generality, NURBS will be considered 
for all future examples. 

Isogeometric analysis relies on the use of the basis functions outlined in this 
section and is conceptually very simple - we use the basis functions used to 
describe the geometry of the problem to approximate the unknown fields in the 
governing PDEs. That is, in the case of elastostatic analysis using BEM, the 
displacement and traction components are approximated using NURBS basis 
functions. The benefit of this approach is clear, since the task of producing 
a boundary discretisation (mesh) is completely provided by CAD and, once 
boundary conditions and material properties have been defined, analysis can be 
carried out immediately. 



3. Conventional BEM 

To illustrate the differences between conventional and isogeometric BEM 
implementations some details of standard BEM technology are presented here. 
The purpose of this section is not to give a complete derivation of the method but 
rather an overview to aid understanding in later sections. The interested reader 
is advised to consult standard BEM references for more details 24 1- 



In this section we describe the direct collocation form of the BEM using piece- 
wise polynomial shape functions; the indirect form of the BEM and the Galerkin 
BEM are not described, though there appears to be no reason why the Galerkin 
BEM cannot be developed in an isogeometric framework. 

To begin, we define the domain of the problem VL with boundary F — dfl. 
We also define two points x' and x, commonly referred to as the source point 
and field point respectively, these points being separated by a distance r given 
by the Euclidean norm 

r:=||x'-x|| (6) 
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(see Fig. The point x' is often referred to as the cohocation point, since in 
the conventional cohocation BEM implementation, the system of equations is 
constructed by taking the collocation point to lie at each nodal point in turn. 
The field point x represents any sampling point, considered in a numerical inte- 
gration scheme, on the portion of boundary over which integration is performed. 
Making the assumption of linear elasticity and in the absence of body forces, 
we can write the displacement boundary integral equation (DBIE) which relates 
displacements and tractions around the boundary F, 

a, (x')^, (x') + £ T,, (x', x)w, (x) dF(x) 

= ^[/,,(x',x)i,(x)dF(x) i,j^l,2 (7) 

where is a jump term that arises from the limiting process of the boundary 
integral on the left hand side of ([7]) and is dependent on the geometry at the 
source point, Uj and tj are the components of displacement and traction around 
the boundary and Uij and are displacement and traction fundamental so- 
lutions relating to a source point direction component i and field point compo- 
nent j. These fundamental solutions for 2D linear elasticity may be found in 
Appendix [Appendix B.l 

In its current form, Eq. ([T]) is not amenable for numerical implementation 
since Uj and tj represent unknown continuous fields. We therefore proceed in 
the usual fashion of discretisation by splitting the boundary of the problem into 
elements with local coordinate rj G [—1,1] over which the geometry and fields 
can be approximated as 

Xe(ry) =^iVf,(77)xb (8) 

b=l 

nt, 

Ue(77) =^iVf,(77H (9) 

b=l 

nt, 

t,(ry) =^iVf,Wt6 (10) 

b=l 

where rih is the number of local basis functions (eg. 6 — 3 for a quadratic 
element), iVf,(?7) are the set of polynomial basis functions (see Fig. [7] for the 
commonly used quadratic Lagrangian basis functions) and Xf,, u;,, t;, are vectors 
of nodal coordinates, displacements and tractions respectively. The subscript e 
has been used in Eqns ([5]) to ([T(I| to denote that the vectors apply to a specific 
element e. By inserting Eqns ^ and (jlOp into (O, the discretised DBIE can be 
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Figure 7: Continuous quadratic basis functions 



written as 



C,,(x')u,(x') 



eb 



e=l b=l 



Nc rib r „ + l 



/ u.A^'Mv))Nbiv)riv)dv 



e^l b^l 




(11) 



where J'^^r]) represents the Jacobian of transformation for element e that maps 
77 — 7> r and 1 < e < A^e is the set of element numbers. 

As mentioned previously, the system of equations is formed by considering 
the collocation point x' to lie at each nodal point in turn. In this way, a set 
of matrices are assembled relating all displacement components and traction 
components. This set of equations can be written as 



with the square matrix H containing all integrals of the kernel plus the jump 
terms Cij, the rectangular matrix G containing all integrals of the Uij kernel, 
and the vectors u, t containing nodal displacement and traction components 
respectively. By prescribing a suitable set of boundary conditions, which may 
consist of a set of displacements and tractions, equation ([T^ may be rearranged 
to give a set of linear equations in the form 



where x represents the vector of unknown degrees of freedom. The linear system 
(fT3|) can be solved for x using conventional or fast solvers 0, while noting 
that A is a full, non-symmetric matrix. 



Hu = Gt 



(12) 



Ax = b. 



(13) 
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Figure 8: Control points and NURBS curve definition for reactor problem. Ap- 
pendix [Xppendi^^ details the control point coordinates and weights. 



4. Isogeometric BEM 

Our attention now focuses on the main idea of the paper: presenting the iso- 
geometric boundary element method implementation for 2D elastostatic analy- 
sis. However, before more details are given, a comment should be made on the 
key difference of IGABEM over conventional BEM. Essentially, it can be reduced 
to the use of NURBS basis functions in place of the conventional polynomial 
counterparts. There are certain consequences for implementation when NURBS 
basis functions are used (such as dealing with nodal points which no longer lie 
on the boundary) , but this key concept should be kept in mind throughout the 
following sections. 

We begin by considering a simple 2D geometry of a nuclear reactor vessel; 
symmetry has been exploited to simplify the problem. The exact geometry can 
be described using the NURBS as illustrated in Fig. [8] All information necessary 
to define the NURBS curve is provided by a CAD model. In this example we 
have chosen, for the entire boundary, to use quadratic basis functions (p = 2) 
which are illustrated in Fig. [S] . 



4.I. Element definition 

As shown in Sec. [31 the boundary integrals in the DBIE are evaluated over 
the entire domain by summing all elemental contributions. But for IGABEM, it 
is not immediately obvious what our definition of elements should be. However, 
we should bear in mind that some notion of 'elements' is used in IGABEM 
simply as a construct for numerical integration, and this is the definition used 
in the remainder of this paper. The element domain is required to cover only 
the portion of the boundary where the relevant basis functions are non-zero. For 
convenience this leads us to a definition of element boundaries as the unique 
values of the knot vector, which can be implemented simply as 
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uniqueKnots = unique (knotVec) ; 

elRanges = [uniqueKnots (1 : end- 1) ' uniqueKnots (2 : end) '] ; 

For example, the first element for the reactor problem is defined by e [0, 1] 
with the set of non-zero basis functions over this element illustrated in Fig. 1101 
For the purposes of numerical integration using Gauss-Legendre quadra- 
ture, local coordinates in the range [—1, 1] over a particular element are most 
commonly used, so a transformation that maps from the parameter space ^ € 
[?iiCa+p+i] to a parent coordinate space defined over an element as ^ e [—1, 1] 
is required. To achieve this, a Jacobian of transformation is used which com- 
prises two terms combined in the chain rule sense: a mapping from the physical 
coordinate space to parameter space (dT/d^) and a mapping from parameter 
space to the local parent coordinate space (d^/d^). This results in the following 
Jacobian of transformation: 

Explicit expressions for the derivatives on the right hand side of Eq. (|14p are 
given in Appendix [Appendix B.2.2| 

4- 1.1. Element connectivity 

Before approximations of the geometry and unknown fields can be given, 
the non-zero basis functions must be determined for a particular element, thus 
forming a connectivity function. If this is done, then a set of local basis functions 
that are related to the global basis functions can be defined as 

N^ii) EE Ra,pm)) (15) 

where the local basis function number 6, element number e and global basis 
function number are related by 

a — conn{e, b) (16) 

where conn{) is a connectivity function. The connectivity function for the reac- 
tor problem is given in Appendix [Appendix C.2[ Using this definition of local 
basis functions, it is now possible to state isogeometric approximations for the 
geometry, displacement and traction as follows: 

P-i-i 

Mi)^Y.^bii)^b (17) 
&=i 
p-i-i 



u,(0 = ^K(Odfc (18) 

&=i 
P+i 

teii)^Y.^biiH (19) 



6=1 
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where Xf,, d;, and qb are vectors of the geometric coordinates, displacement 
coefficients and traction coefficients respectively, associated with the control 
point corresponding to the basis function b. It should be noted that we have 
used the term coefficient since the NURBS basis functions do not necessarily 
obey the Kronecker-delta property (in contrast to the polynomial basis functions 
shown in Fig. [7]) . Therefore, the terms and do not necessarily represent 
real displacements and tractions. 

Is it is now a simple case of substituting Eqns (|18p and ([T5| into the DBIE 
of ([7]) while using the Jacobian given by ([Tl)). This results in the following 
discretised equation for IGABEM: 



1=1 



EE 

e=l 1=1 

We P+1 

EE 

e=l 1=1 



T,,(x',x(0)7Vf(0J(ad| 
C/,,(x',x(|))iVf(OJ(Od| 



1^ 



(20) 



where d''f 



and represent components of the vectors db and qt for the element 



e. We denote e as the element containing the collocation point x' and C is the 
local coordinate of the collocation point in element e. The reason for specifiying 
these terms concerns the term Uj(x') as seen in Eq. ([TT]) . In the conventional 
implementation where collocation occurs at nodal points, the Kronecker-delta 
property of the basis functions ensures that at the collocation point itself the 
basis functions are interpolatory. In contrast to this, the IGABEM formulation 
cannot guarantee that this is true and the displacement must be interpolated 
as ^^,(x') = Et>f(l')rf^ 

4-2. Collocation point definition 

A significant change in IGABEM over conventional BEM is in the location of 
collocation points, since the normal practice of collocation at nodal positions is 
no longer valid. We can easily see why this is the case by inspecting the position 
of the control points in Fig. [5] it is evident that there is one control point in this 
example that does not lie on the boundary. Indeed, for most curved boundaries 
the control points will not lie on the boundary. Since the control points can be 
interpreted as nodes in the IGABEM formulation, this presents a problem since 
it is essential that collocation takes place on the boundary F. To overcome this, 
we choose to use the Greville abscissae definition 17| , [20[ to define the position 
of collocation points in parameter space. This is defined as: 



A/ Ca+l + Ca+2 + • • • + Ca+P , o 



(21) 



In Matlab@ , this is easily evaluated as: 
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Figure 11: Collocation point and element definitions for reactor problem. 



collocPts = zeros(l,n); 
for i=l:n 

collocPts(i) = sum(knotVec(i+l : i+p) ) / p; 

end. 

where n denotes the number of control points, p is the curve order and knotVec 
is the knot vector. If this definition is applied to the geometry of the reactor 
vessel, then the collocation points are as shown in Fig. [11] The coordinates in 
physical space can be found by using (HI) with the NURBS basis functions given 
by (O. The definition of the element boundaries is also illustrated in this figure. 

4.3. Implementation aspects 

The previous two sections outlining the definition of boundary elements and 
collocation points within an isogeometric BEM framework specify the key as- 
pects that must be changed for an existing collocation BEM implementation. 
Our attention now focuses on the implementation for the reactor vessel problem 
to illustrate these aspects. 



4.3.1. Refinement 

The discretisation in Fig. |5| represents the geometry of the reactor prob- 
lem exactly, but experienced analysts will be acutely aware that, although the 
geometry may be captured, the basis functions may be insufficient to capture 
the gradients in the unknown fields, consequently leading to large errors in the 
solution. This leads to one of the most beneficial properties of B-splines and 
NURBS for numerical analysis, which is that the mesh can be refined to ar- 
rive at a richer set of basis functions while preserving the exact geometry at all 
stages. In the current paper we present two types of refinement: knot insertion 
and order elevation, termed h-refinement and p-refinement in numerical analysis 
literature. The algorithms for these refinement processes are standard in CAD, 
with details given in 24|, 26| and source code given in [l|. 
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Figure 12: Control points and NURBS curve definition for reactor problem after 
knot insertion. 




Figure 13: NURBS basis functions for reactor problem after uniform knot insertion 
(h-refinment) 



We can illustrate both knot insertion and order elevation using the reactor 
problem shovi^n in Fig. [8] Figs [12] and |T3| illustrate knot insertion where addi- 
tional knots have been inserted uniformly into the original knot vector. Figs[T4l 
and [Tsj illustrate order elevation where the basis functions have been increased 
from quadratic {p — 2) to cubic (p = 3). Both types of refinement introduce 
changes to the knot vector and set of control points. A third type of refine- 
ment which is often referred to as k-refinement can also be considered which is 
a combination of p- and h- refinement. It arises due to the non-commutative 
nature of the aforementioned refinements, k-refinement is not considered in the 
present paper. 
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4-3.2. Integration 

A key feature of any BEM implementation is the evaluation of the boundary 
integrals containing the kernels over element domains. It is well-known that 
both regular and singular integrands are found depending on the position of 
the collocation point relative to the field element. Essentially, the evaluation of 
BEM integrals is split into three different types described as 

1. Regular integration: the collocation point lies in an element different 
from the field element. 

2. Nearly singular integration: the collocation point lies in a element not 
on but near the field clement. 

3. Singular integration: the collocation point lies in the field element and 
can be one of two types: 

Strongly singular integrals: (T!y kernel, 0{l/r) in 2D) 
Weakly singular integrals: {Uij kernel, 0(ln(l/r)) in 2D) 

In the present study, we choose to treat the regular and nearly-singular integrals 
in the same manner, although several methods exists for the efficient treatment 
of integrals of the latter type [32,],^2l[, 16 1. The evaluation of singular integrals 



must be given close consideration since they are found to have a large infiuence 
on the accuracy of the resulting solution. 

The present work uses the subtraction of singularity method (SST) [3] to 
evaluate strongly singular integrals, in which the integrand is split into its reg- 
ular and singular parts. Further details of the method are shown in [Soj with 
a full implementation given in [ij]. The idea is to separate the integral into a 
regular part (which can evaluated using standard quadrature routines) and a 
singular part which can be evaluated analytically (for 2D boundary integrals). 
If the SST method is used, however, the jum p te rm C'ij must be calculated ex- 
plicitly (c.f. the rigid body motion technique [l5| which calculates it implicitly). 
However, a simple formula exists [3] for the 2D case which is restated here as 



C ' 



87r(l - P) 



4{1 - D)i9i - 62) 
+ (sin26li - sin 26*2) 

cos 26*2 — cos 201 



cos 202 — cos 29i 

4(1 - 9){0, ~ 62) 
- (sin26li - sin 26*2) 



(22) 



where v = v for plane strain and v = v ji^ ^ v) for plane stress and the angles 
0i,^are related to the normals of the surface at the collocation point (see Fig. 2 
in Q). 

For weakly singular integrals, a variety of techniques are available includ- 
ing specific logarithmic quadrature points and weights (30| . but in the present 
study we choose to use the Telles transformation [s^] which cancels the sin- 
gularity leaving a regular integrand. This is achieved through the following 
transformation: 

, +Y(y^+3) 

^ (l + 37'2) ^^^^ 
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where 



i = V^'(^" - 1) + 1^" - 1| + V - 1) - 1^" - 1| + ^' (24) 

^' denotes the location of the smgularity in the parent space S [—1,1]) and 7 
represents the new integration variable. Therefore, a Jacobian which transforms 
from the parent space ^ G 1] to 7 is required. This is given as: 



■J«=4^l7 (25) 



Using Eqn. (P^. (1^^ and (l?5]) . the transformation of the integration can be 
expressed as: 



-1 f+i 

/m- / / 
1 J -I 



(7-Yr + y(y^ + 3) 

(1 + 37'2) 



3(7-7') 



l\2 



1 + 37 



7fd7 (26) 



A derivation of the transformation is given completely in [sS] with a full imple- 
mentation for the Uij kernels given in [1]. 

i.3.3. IGABEM algorithm 

Now that the integration routines which form the core of IGABEM imple- 
mentation have been outlined, we are in a position to give an overview of the 
entire IGABEM algorithm, illustrated in Algorithm [T] 

In this algorithm the H and G matrices have been calculated explicitly before 
boundary conditions are applied to arrive at the final system of equations. For 
efficiency, commercial BEM implementations often calculate A and z directly, 
thereby making the construction of the H and G matrices redundant. 

^.4. Example 

Finally, an example problem using the geometry of the reactor used through- 
out the current paper is defined and analysed. The problem geometry, bound- 
ary conditions and material properties are shown in Fig. 1161 where symmetrical 
boundary conditions are applied at cc = 100 and y = 0, and a constant pressure 
is exerted on the inner boundary. Plane strain is assumed. 

By applying knot insertion to the original mesh shown in Fig. Illl the dis- 
cretisation shown in Fig. [T7] was used to perform the IGABEM analysis. The 
results are shown in Fig. [18] by plotting the exaggerated displacement profile 
with FEM results also plotted for comparison. As can be seen, excellent agree- 
ment is obtained. In addition, a convergence study was carried out to assess the 
accuracy of IGABEM over a standard BEM implementation with quadratic ba- 
sis functions. Using h-refinement in both methods, the L2 norm in displacement 
was calculated for each mesh as: 

M\L.^^ I EWdr (27) 
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Algorithm 1 IGABEM algorithm 



Read CAD input data > e.g. Control points, knot vector 

Read material properties, boundary conditions 

Perform mesh refinement > e.g. knot insertion or order elevation 



for c ^ 1, rir do 



> Loop over collocation pts 



jumpTerm ^ calcJ umpT erm{collocN ormals{c)) \> Calculate Ci- 



for e l,nei do 

elementConn ^ glhConn{e) 
elRange ^ ElmtRanges{e) 



Loop over elements 
> Element connectivity 
> Range of element 



if c e elRange then 

Hs„{, SSTIntegration{c,e,Cij) 
Gsub ^ TellesIntegration{c, e) 

else 

[Hs„b, Gs„b] -e- GaussLegendreQuad{c,e) 
end if 



> Singular integration 
> SST integration 
> Telles transformation 
Non-singular integration 



H(c, elementConn) i— Hsub > Add submatrix to global H matrix 
G(c, elementConn) Ggufe > Add submatrix to global G matrix 
end for 
end for 



[A,z] <— apply BoundConds{H,G) 

X ^ sol (Y=(A. z) 



> Apply boundary conditions 

[> Sohc s\ st.('m of (Xjuations 
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Figure 16: Definition of nuclear reactor geometry, boundary conditions and material 
properties. 

The results obtained are shown in Fig. [19] where a significant improvement in 
accuracy over the standard BEM implementation can be seen. The reference 
solution corresponds to the converged result of a BEM analysis using quadratic 
basis functions. 

5. Conclusions 

The implementation aspects of an isogeometric BEM were outlined, with 
attention paid to areas which differ from a conventional BEM implementation. 
What is evident is that IGABEM present a particularly attractive approach for 
analysis, since the data provided by CAD can be used directly without the need 
to create a mesh. In addition, refinement schemes were outlined which provide 
more refinement or richer basis functions in required areas. Finally, an example 
was used to illustrate the simplicity and accuracy of the method where it was 
shown that good aggrement with FEM was obtained, and in addition, significant 
improvements in accuracy over a standard BEM implementation with quadratic 
basis functions were demonstrated. 
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Figure 17: Mesh used for reactor problem IGABEM analysis detailing element edges 
and collocation point positions. 




Figure 18: 



Comparison of IGABEM and FEM results for reactor problem - exagger- 
ated displacement profile. 



22 



.418 - 




.408 - [. 



1000 2000 3000 4000 

DOF 

Figure 19: Comparison of L2 displacement norms for reactor problem using IGA- 
BEM and standard BEM with quadratic basis functions. 



Appendix A. B-splines/NURBS 

Appendix A.l. B- spline derivatives 

The first order derivative of the B-spline basis function is expressed as 

;^iV,,p(e) - ■ ^ ■ A^a.p-i(C) - 7 ^^A^a+i,p-i(0- (A.l) 

where p is the polynomial order, a the basis function index. The higher order 
derivatives can be obtained by differentiating the two sides of equation (jA.ip : 

(A.2) 

From 

Eqn(|XT])and(IX2l), 

we can express high order derivatives with Na^p-k, • • • , Na+k,p-k'- 
^^-^■^^^^ = (^^Eafc.;>iVa+..p-.(0, (A.3) 



with 



ao,o = 1, 



ak,o 



ak,b = -z 6=l,...,fc~l, 

t,a+p+b-k+l — t,a+b 



ak,k 



— Qfc-l,fc-l 
ia+p+l — S,a+ 
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Appendix A. 2. NURBS derivatives 

From Eqn ([S]), we can give the first order derivative of NURBS basis function 

where = f^N^^O and 

n 

W^'(e)=E^a,p(eK. (A.5) 

a=l 

We introduce some notations for convenience 

d^ 

^i'^iO^Wa^NaAO, (no sum on z) (A.6) 

and 

W^'HO^^WiO- (A.7) 

Then higher-order derivatives of these rational functions may be expressed in 
terms of lower-order derivatives as 

-^RaiO = , (A.8) 



where 

'k\ k\ 



hj b\{k-h)\ 



(A.9) 



Appendix B. BEM 



Appendix B.l. 2D Fundamental solutions 

Denoting ^ as the shear modulus, v as Poisson's ratio, 5ij as the kronecker- 
delta function defined by 

%4 ] (B.l) 

and noting that a comma denotes differentiation, the fundamental solutions for 
2D linear elasticity are given by: 

t/,(x',x) . {(3-4.)ln (i) (B.2) 

-1 f dr 

T„-(x',x) = 4^(-^ _ \d^^^^ ^^^^^^ + 2r ^'"j] - (1 - 2zy)(r - r^n,;) 

(B.3) 
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Appendix B.2. Boundary element integration parameters 
Appendix B.2.1. Normals 

The normals can be calculated by: 



1 



with 



rix = 
ny = 

MO 

MO 



JiO 
1 



dyiO 
MO 



dC 



dC 

d^ d^ 
Appendix B.2.2. Jacobian of transformation 

dm _o-^s 



b=l 
p+1 

E 

6=1 



-Vb- 



(B.4) 
(B.5) 

(B.6) 
(B.7) 



2 

where and denotes the values of knots at the beginning and end 
element respectively (assuming an outward pointing normal). 



dr 

d? 



(B.8) 
of the 

(B.9) 
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Appendix C. Reactor problem data 



Appendlr C.I. C'oiil rol poiiil.'i (vitd wc/ifjlil.'i 



index (a) 


Control point coordinate (Pa) 


Weight Wa 


1 


(0, 0) 


1 


2 


(20, 0) 


1 


3 


(40, 0) 


1 


4 


(40, 60) 


V2/2 


5 


(100, 60) 


1 


6 


(100, 80) 


1 


7 


(100, 100) 


1 


8 


(72.5, 100) 


1 


9 


(45, 100) 


1 


10 


(45, 87.5) 


1 


11 


(45, 75) 


1 


12 


(35, 75) 


1 


13 


(25, 75) 


1 


14 


(25. 57.5) 


1 


15 


(25. 40) 


1 


16 


(17.5, 40) 


1 


17 


(10, 40) 


1 


18 


(10, 27.5) 


1 


19 


(10. 15) 


1 


20 


(5, 15) 


1 


21 


(0, 15) 


1 


22 


(0, 7.5) 


1 


23 


(0, 0) 


1 



Appendix G.2. Connectivity matrix 



element index (e) 


global basis index (a) for (61,62,63) 


1 


1, 2, 3 


2 


3, 4, 5 


3 


5, 6,7 


4 


7, 8,9 


5 


9, 10, 11 


6 


11, 12, 13 


7 


13, 14, 15 


8 


15, 16, 17 


9 


17, 18, 19 


10 


19, 20, 21 


11 


21, 22, 1 



V 
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